1. Conditional Mean and Variance

Suppose $Y_t$ is AR(1) series observation
(unconditional) Mean \[ E( Y_t ) = 0 \]
(unconditional) Variance \[ V( Y_t ) = \gamma(0) = (1+\phi_1^2) \sigma^2 \]

Conditional Mean

Conditonal Mean:
\[ E\Big(Y_t \hspace2mm \Big| \hspace2mm \mbox{all variables realized by yesterday. }\Big) \]
Conditonal mean of AR(1): $Y_t = Y_{t-1} + e_t $ \[ \begin{align} E\Big(Y_t \hspace2mm \Big| \hspace2mm Y_{t-1}, Y_{t-2}, \ldots, e_{t-1}, e_{t-2}, \ldots\Big) &= E \Big(\phi Y_{t-1} + e_t \hspace2mm \Big| \hspace2mm Y_{t-1}, Y_{t-2}, \ldots, e_{t-1}, e_{t-2}, \ldots\Big)\\ \\ &= \phi Y_{t-1} + E(e_t) \hspace2mm = \hspace2mm \ph Y_{t-1}. \end{align} \]
Conditonal variance \[ \begin{align} V\Big(Y_t \hspace2mm \Big| \hspace2mm Y_{t-1}, Y_{t-2}, \ldots, e_{t-1}, e_{t-2}, \ldots\Big) &= V\Big(\phi Y_{t-1} + e_t \hspace2mm \Big| \hspace2mm Y_{t-1}, Y_{t-2}, \ldots, e_{t-1}, e_{t-2}, \ldots\Big)\\ \\ &= V(e_t) \hspace2mm = \hspace2mm \sigma^2 \end{align} \]

    n=200
    Y <- arima.sim(n = n, list(ar = c(0.8) ))
    
    v.AR1 = (1+.8^2)*1
    
    plot(Y, type="o", main="AR(1)", xlim=c(0,n*1.1))
    abline(h=0)
    abline(h=c(1,-1)*1.96*v.AR1, col="red")

    lines(n+1, .8*Y[n], type="p", col="blue")
    lines(c(n+1,n+1), .8*Y[n]+c(1.96,-1.96), type="p", col="red")

Suppose $Y_t$ is MA(1) series observation
(unconditional) Mean \[ E( Y_t ) = 0 \]
(unconditional) Variance \[ V( Y_t ) = \gamma(0) = (1+\theta_1^2) \sigma^2 \]
Conditonal mean of MA(1): $Y_t = e_t + \theta_1 e_{t-1}$ \[ \begin{align} E\Big(Y_t \hspace2mm \Big| \hspace2mm Y_{t-1}, Y_{t-2}, \ldots, e_{t-1}, e_{t-2}, \ldots\Big) &= E \Big( e_t + \theta_1 e_{t-1} \hspace2mm \Big| \hspace2mm Y_{t-1}, Y_{t-2}, \ldots, e_{t-1}, e_{t-2}, \ldots\Big) \\\\ &= E ( e_t )+ \theta_1 e_{t-1} \hspace2mm = \hspace2mm \theta_1 e_{t-1} \end{align} \] Note that $e_{t-1}$ is not observable.
Conditonal variance of MA(1): $Y_t = e_t + \theta_1 e_{t-1}$ \[ \begin{align} Var\Big(Y_t \hspace2mm \Big| \hspace2mm Y_{t-1}, Y_{t-2}, \ldots, e_{t-1}, e_{t-2}, \ldots\Big) \ &= Var\Big( e_t + \theta_1 e_{t-1} \hspace2mm \Big| \hspace2mm Y_{t-1}, Y_{t-2}, \ldots, e_{t-1}, e_{t-2}, \ldots\Big)\\ &= Var( e_t ) \hspace2mm = \hspace2mm \sigma^2. \end{align} \]
AR(1) \[ \mbox{ Uncond'l } \hspace20mm \mbox{ cond'l } \\ E(Y_t) = 0, \hspace28mm E(Y_t|\omega_{t-1}) = \phi_1 Y_{t-1}, \\ Var(Y_t) = (1+\phi_1^2)\sigma^2 \hs10mm Var(Y_t|\omega_{t-1}) = \sigma^2 \]
MA(1) \[ \mbox{ Uncond'l } \hspace20mm \mbox{ cond'l } \\ E(Y_t) = 0, \hspace28mm E(Y_t|\omega_{t-1}) = \theta_1 e_{t-1}, \\ Var(Y_t) = (1+\th_1^2)\sigma^2 \hs10mm Var(Y_t|\omega_{t-1}) = \sigma^2 \]
For ARMA(p,q) model, conditional mean changes, but conditional variance is constant.

  x  <- arima.sim(n = 50, list(ar = c(0.8) )) 
  plot(x, type="o", main="AR(1)"); abline(h=0)

  x  <- arima.sim(n = 50, list(ma = c(0.5) )) 
  plot(x, type="o", main="MA(1)"); abline(h=0)

Financial Return

\[ \begin{align} X_t &= \mbox{ Stock Price (observation) } \\ \\ Y_t &= \ln(X_t) - \ln(X_{t-1}) \hs5mm : \mbox{ log return } \end{align} \]

  library(quantmod)
  source("http://gozips.uakron.edu/~nmimoto/477/TS_R-90.txt")

  getSymbols("AAPL")    #- download from Yahoo!

## [1] "AAPL"

  Apple <- Ad(AAPL)['2007::2009']

  plot(Apple)

  plot(log(Apple))

  plot(diff(log(Apple)))

  Randomness.tests(as.numeric(diff(log(Apple))))

##   B-L test H0: the sereis is uncorrelated
##   M-L test H0: the square of the sereis is uncorrelated
##   J-B test H0: the sereis came from Normal distribution
##   SD         : Standard Deviation of the series

##      BL15  BL20  BL25 ML15 ML20 JB    SD
## [1,] 0.07 0.006 0.008    0    0  0 0.028

  getSymbols("^GSPC")    #- SP500 download from Yahoo!

## [1] "GSPC"

  SP500 <- Ad(GSPC)['2017::']

  plot(SP500, type="o")

  plot(log(SP500), type="o")

  plot(diff(log(SP500)), type="o")

  Randomness.tests(as.numeric(diff(log(SP500))))

##   B-L test H0: the sereis is uncorrelated
##   M-L test H0: the square of the sereis is uncorrelated
##   J-B test H0: the sereis came from Normal distribution
##   SD         : Standard Deviation of the series

##       BL15  BL20  BL25 ML15 ML20 JB    SD
## [1,] 0.001 0.001 0.003    0    0  0 0.007

Stylized Facts about Financial Return

Uncorrelated
Squares are correlated
Clustering
Asymmetry
Heavy Tailed unconditional and conditional distribution

Heteroschedasticity

Don’t confuse the conditional heteroscedasticity with (unconditonal) heteroscedasticity:

  D <- read.csv("http://gozips.uakron.edu/~nmimoto/pages/datasets/gas.csv")
  D1 <- ts(D[,2], start=c(1956, 1), freq=12)
  plot(D1, type='o')

2. ARCH model

Engle (1985) AutoRegressive Conditionally Heteroscedastic Model \[ \begin{align} Y_t &= \sigma_t e_t \hs10mm e_t \sim_{iid} \cN(0,1) \\\\ \sigma_t^2 &= \omega+ \alpha Y_{t-1}^2 \end{align} \]
Mean \[ E(Y_t) = \sigma_t E(e_t) = 0 \]
Variance \[ V(Y_t) = V(\sigma_t)V(e_t) = \]

Conditonal Mean and Variance of ARCH

Conditional mean \[ E[Y_t \Big| Y_{t-1}, e_{t-1}, \ldots] \hspace3mm = \hspace3mm \sigma_t E[e_t] = 0 \]
Conditional variance \[ V[Y_t \Big| Y_{t-1}, e_{t-1}, \ldots] \hspace3mm = \hspace3mm \sigma_t^2 V[e_t] = \sigma_t^2 \]

ARCH is uncorrelated

$Y_t$ is uncorrelated, so it will pass the Ljung-Box test.
But $Y_t^2$ is correlated, so it will NOT pass the McLeod-Li test.

3. GARCH model

GARCH(1,1) model \[ \begin{align} Y_t &= \sigma_t e_t \hspace10mm e_t \sim_{iid} N(0,1) \\\\ \sigma_t^2 &= \omega + \alpha Y_{t-1}^2 + \beta \sigma^2_{t-1} \end{align} \]

Conditional Mean and Var of GARCH

Conditional Mean: $0$
Conditional Variance : $\sigma_t^2$

Example: Daily SPY

Daily Price of SP500 ETF (SPY) from Jan 02 2000 to Dec 31 2014

  D <- read.csv("http://gozips.uakron.edu/~nmimoto/pages/datasets/SPY.csv", header=T)
  source("http://gozips.uakron.edu/~nmimoto/477/TS_R-90.txt")
    
  head(D)

##   ixDay      Date   Weekday X SPY.Open SPY.High SPY.Low SPY.Close
## 1     1  1/3/2000    Monday 1   148.25   148.25  143.88    145.44
## 2     2  1/4/2000   Tuesday 2   143.53   144.06  139.64    139.75
## 3     3  1/5/2000 Wednesday 3   139.94   141.53  137.25    140.00
## 4     4  1/6/2000  Thursday 4   139.62   141.50  137.75    137.75
## 5     5  1/7/2000    Friday 5   140.31   145.75  140.06    145.75
## 6     6 1/10/2000    Monday 6   146.25   146.91  145.03    146.25
##   SPY.Volume SPY.Adjusted
## 1    8164300       110.33
## 2    8089800       106.01
## 3   12177900       106.20
## 4    6227200       104.50
## 5    8066500       110.57
## 6    5741700       110.94

  X <- ts(D[,7], start=1)
  plot(X)

  X1 <- diff(log(X))
  plot(X1)

  layout(matrix(1:2, 1, 2))
  acf(X1)
  acf(X1^2)

  layout(1)


  # Fit GARCH model 
  library(fGarch)
  Fit1 <-  garchFit(~ garch(1,1), data=X1, cond.dist="norm",  include.mean = FALSE, trace = FALSE)
    
  ##  Fit1@fit$par          #-- estimated parameters
  ##  Fit1@residuals        #-- this is not the garch residuals!!!! same as X1
  ##  Fit1@sigma.t          #-- estimated sig_t
  ##  Fit1@residuals/Fit1@sigma.t       #-- this is the (standardized) GARCH residuals

  print(Fit1@fit$ics)   #-- AIC and BIC are here

##       AIC       BIC       SIC      HQIC 
## -6.470319 -6.465360 -6.470321 -6.468556

  plot(X1)

  plot(Fit1@sigma.t, type="l")   #- estimated Sigma_t

  plot(X1)
  lines( 1.96*Fit1@sigma.t, type="l", col="red")
  lines(-1.96*Fit1@sigma.t, type="l", col="red")

  res1  <- X1/Fit1@sigma.t   #-- this is the residuals
  Randomness.tests(res1)

##   B-L test H0: the sereis is uncorrelated
##   M-L test H0: the square of the sereis is uncorrelated
##   J-B test H0: the sereis came from Normal distribution
##   SD         : Standard Deviation of the series

##      BL15 BL20 BL25  ML15  ML20 JB    SD
## [1,]    0    0    0 0.719 0.855  0 1.001

Representing $\sigma^2_t$ using $Y_t^2$

\[ \begin{align} \sigma_t^2 &= \omega + \alpha Y_{t-1}^2 + \beta \sigma^2_{t-1} \\\\ &= \omega + \alpha Y_{t-1}^2 + \beta \Big(\omega + \alpha Y_{t-2}^2 + \beta \sigma^2_{t-2} \Big)\\ \\ &= \omega + \beta \omega + \alpha Y_{t-1}^2 + \beta \alpha Y_{t-2}^2 + \beta^2 \sigma^2_{t-2} \\ \\ &= \omega + \beta \omega + \alpha Y_{t-1}^2 + \beta \alpha Y_{t-2}^2 + \beta^2 \Big(\omega + \alpha Y_{t-3}^2 + \beta \sigma^2_{t-3} \Big) \\ \\ &= \omega + \beta \omega + \beta^2 \omega + \alpha Y_{t-1}^2 + \beta \alpha Y_{t-2}^2 + \beta^2 \alpha Y_{t-3}^2 + \beta^3 \sigma^2_{t-3} \end{align} \]

Continuing, we get \[ \begin{align} &= \omega (1 + \beta + \beta^2 + \cdots ) + \alpha \sum_{i=0}^k \beta^i Y_{t-1-i}^2 + \beta^{k+1} \sigma^2_{t-1-k} \\\\ &= \frac{\omega}{1-\beta} + \alpha \sum_{i=0}^\infty \beta^i Y_{t-1-i}^2 \end{align} \]

\[ \sigma_t^2 \hspace3mm = \hspace3mm \frac{\omega}{1-\beta} + \alpha \sum_{i=0}^\infty \beta^i Y_{t-1-i}^2 \] We will use the truncated and estimated version of this, \[ \hat \sigma_t^2 \hspace3mm = \hspace3mm \frac{\hat \omega}{1- \hat \beta} + \hat \alpha \sum_{i=0}^{t-1} \hat \beta^i Y_{t-1-i}^2 \]

Residuals of GARCH process

\[ Y_t = \sigma_t e_t \]

Using observation $\{Y_1, \ldots, Y_n\}$, the residuals are \[ \hat e_t = Y_t/ \hat \sigma_t \]

Quasi-Maximum Likelihood Estimator

Importatnce of Conditional Distribution

\[ \begin{align} Y_t &= \sigma_t e_t \\\\ \sigma_t^2 &= \omega + \alpha Y_{t-1}^2 + \beta \sigma^2_{t-1} \end{align} \]

[Conditonal distribution] $=$ [Distribution of $e_t$]
Guessing correct distribution for $e_t$ is very important in GARCH parameter estimation.

Built-in distributions

  x <- seq(-10,10, .1)
    
  #- Standardized t-distribution 
  plot(x, dstd(x, mean = 0, sd = 1, nu = 5), type="l" )
  lines(x, dnorm(x, mean = 0, sd = 1), col="red" )

  #- Generalized Error Distribution
  plot(x, dged(x, mean = 0, sd = 1, nu = 4), type="l" )
  lines(x, dnorm(x, mean = 0, sd=1) , col="red" )

  #- Skewed-Standardized t-distribution 
  plot( x, dsstd(x, mean = 0, sd = 1, nu = 5, xi = 1.5), type="l" )
  lines(x, dnorm(x, mean = 0, sd = 1), col="red" )

  #- Skewed-Generalized Error Distribution
  plot(x, dsged(x, mean = 0, sd = 1, nu = 5, xi = 1.5) , type="l" )
  lines(x, dnorm(x, mean = 0, sd = 1), col="red" )

#--- 5. Simulation and estimation of GARCH -------------------------------------------------
  n     <- 1000
  theta <- c(.024, .1, .8)   #- Parametes of GARCH(1,1)


  #- Specify parameters and conditional distribuiton (distribution of e_t)
  spec <- garchSpec(model = list(omega=theta[1], alpha=theta[2], beta=theta[3]),  cond.dist="norm")

  #- Simulate GARCH using above spec
  x    <- garchSim(spec, n = n, extended=FALSE)

  #- Estimate parameters of GARCH 
  Fit1 <-  garchFit(~ garch(1,1), data=x, cond.dist="norm",  include.mean = FALSE, trace = FALSE)
  coef(Fit1) #- estimated parameters


  #--- Monte Carlo Estimation to see parameter MSE --
  n     <- 1000
  theta <- c(.1, .1, .8)   #- Parametes of GARCH(1,1)


  #- Specify parameters and conditional distribuiton (distribution of e_t)
  spec <- garchSpec(model = list(omega=theta[1], alpha=theta[2], beta=theta[3]),  cond.dist="norm")


  Param <- matrix(0, 1000, 3)
  for (i in 1:1000) {

    print(i)

    Y    <- arima.sim(

    x    <- garchSim(spec, n = n, extended=FALSE)
    Fit1 <- garchFit(~ garch(1,1), data=x, cond.dist="norm",  include.mean = FALSE, trace = FALSE)

    Param[i,] <- coef(Fit1) #- estimated parameters

  }
  
  SE <- (Param - t( matrix(theta, 3, 1000) ))^2
  MSE <- apply(SE, 2, mean) #- MSE of each parameter
  rMSE <- sqrt(MSE);  rMSE

Simulation Study

True parameter (.024, .1, .8)

True cond’l dist: Normal Estimated using: Normal
True cond’l dist: std(5) Estimated using: Normal
True cond’l dist: sged(skew=.7, shape=1.45) Estimated using: Normal

  MSE1/MSE1   #- Relative MSE 
  [1] 1 1 1

  MSE2/MSE1
  [1] 1.584734 3.396001 1.831015

  MSE3/MSE1
  [1] 1.198293 1.613116 1.355613

Lecture 25 - GARCH model

Nao Mimoto - Dept. of Statistics : The University of Akron

Apr 6 2018